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ABSTRACT 

Moment  models  of  carrier  transport,  derived  from  the  Boltzmann  equation,  have  made 
possible  the  simulation  of  certain  key  effects  through  such  realistic  assumptions  as  energy 
dependent  mobility  functions.  This  type  of  global  dependence  permits  the  observation  of 
velocity  overshoot  in  the  vicinity  of  device  junctions,  not  discerned  via  classical  drift-diffusion 
models,  which  are  primarily  local  in  nature.  It  has  been  found  that  a  critical  role  is  played 
in  the  hydrodynamic  model  by  the  heat  conduction  term.  When  ignored,  the  overshoot  is 
inappropriately  damped.  When  the  standard  choice  of  the  Wiedemann- Franz  law  is  made  for 
the  conductivity,  spurious  overshoot  is  observed.  Agreement  with  Monte-Carlo  simulation  in 
this  regime  has  required  empirical  modification  of  this  law,  as  observed  by  IBM  researchers, 
or  nonstandard  choices.  In  this  paper,  simulations  of  the  hydrodynamic  model  in  one  and 
two  dimensions,  as  well  as  simulations  of  a  newly  developed  energy  model,  the  RT  model,  will 
be  presented.  The  RT  model,  intermediate  between  the  hydrodynamic  and  drift-diffusion 
model,  was  developed  at  the  University  of  Illinois  to  eliminate  the  parabolic  energy  band  and 
Maxwellian  distribution  assumptions,  and  to  reduce  the  spurious  overshoot  with  physically 
consistent  assumptions.  The  algorithms  employed  for  both  models  are  the  essentially  non- 
oscillatory  shock  capturing  algorithms,  developed  at  UCLA  during  the  last  decade.  Some 
mathematical  results  will  be  presented,  and  contrasted  with  the  highly  developed  state  of 
the  drift-diffusion  model. 
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1  Introduction 

During  the  last  decade,  device  modeling  has  attempted  to  incorporate  general  carrier  heating, 
velocity  overshoot,  and  various  small  device  features  into  carrier  simulation.  The  popular 
wisdom  emerging  from  such  concentrated  study  holds  that  global  dependence  of  critical 
quantities,  such  as  mobilities,  on  energy  and/or  temperature,  is  essential  if  such  phenomena 
are  to  be  modeled  adequately.  In  this  paper,  we  examine  in  detail  the  simulation  of  two 
such  energy  models,  including  the  hydrodynamic  model  and  the  RT  model.  We  describe  the 
models,  summarize  some  associated  mathematical  results,  as  well  as  the  basic  features  of  the 
numerical  algorithm,  and  then  present  the  results  of  extensive  numerical  simulations  for  two- 
dimensional  MESFET  devices,  and  for  one-dimensional  diodes.  Both  models  represent  one 
carrier  flow.  The  hydrodynamic  model  contains  hyperbolic  modes  related  to  the  momentum 
equations,  while  the  RT  model  does  not  possess  such  modes.  In  both  cases,  however,  we 
employ  a  conservation  law  format,  and  numerical  methods  suitable  for  such  systems.  The 
ENO  (essentially  non-oscillatory)  method  employed  makes  use  of  adaptive  stencils,  and  is 
particularly  adept  at  shock  capturing  if  the  parameter  regime  crosses  from  supersonic  to 
subsonic.  Even  if  this  does  not  occur,  the  convective  terms  are  effectively  discretized,  via 
this  procedure,  in  both  models.  The  first  use  of  such  methods  in  device  simulation  was  in 
[7],  followed  by  the  study  [6],  in  which  shocks  were  detected  in  micron  devices  at  liquid  Ni¬ 
trogen  temperatures,  and  at  room  temperature  in  shorter  devices,  by  independent  numerical 
techniques. 

Our  development  of  the  RT  model  follows  that  of  [5].  These  researchers  attempted 
to  utilize  a  microscopic  relaxation  time  approximation,  which  would  allow  for  nonparabolic 
energy  bands  and  non-Maxwellian  distribution  functions.  The  approach  allows  for  parameter 
fitting  of  certain  key  quantities  via  Monte-Carlo  simulation. 

One  of  the  principal  conclusions  of  the  paper  is  the  essential  dependence  of  the  hydrody¬ 
namic  model  upon  the  heat  conduction  term.  Standard  choices  lead  to  numerically  detected 
spurious  overshoot  at  the  drain  junction  of  an  —  n  —  n'^  diode,  while  other  choices  signifi¬ 
cantly  damp  this  overshoot.  Monte-Carlo  simulations  show  that  substantial  underestimation 
occurs  when  the  heat  conduction  term  is  neglected.  We  refer  the  reader  to  [10],  and  to  the 
simulation  results  of  this  paper  for  amplification.  The  RT  model  was  developed,  partly  in 
response  to  the  continuing  debate  concerning  heat  conduction  processes  in  the  hydrodynamic 
model. 

The  status  of  mathematical  results  differs  sharply  between  the  hydrodynamic  model, 
on  the  one  hand,  and  the  drift-diffusion  model  on  the  other.  For  the  former,  we  hav’e 
summarized  two  results,  one  by  Gamba  (cf.  [8])  for  an  idealized  model,  in  which  the  adi- 
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abatic  relation  is  employed,  and  another  by  Gardner,  Jerome,  and  Rose  (cf.  [9])  in  which 
a  Newton- Kantorovich  theorem  is  developed  for  the  n'*'  —  n  —  diode,  yielding  existence 
and  convergence  in  a  specialized  subsonic  regime.  The  drift-diffusion  model,  on  the  other 
hand,  has  been  widely  studied.  Existence  and  approximation  results  have  been  carefully 
developed,  although  uniqueness  is  still  not  well  understood  for  this  model.  Existence  for 
the  steady-state  model  is  due  in  varying  degrees  of  generality  to  many  authors,  including 
Mock  ([17]),  Seidman  ([21]),  and  the  first  author  ([13]).  A  convergence  theory,  based  upon 
a  calculus  due  to  Krasnosel’skii,  was  presented  in  ([15]).  Mathematical  results  have  not  yet 
been  developed  for  the  strongly  nonlinear  RT  model. 


2  Hydrodynamic  and  Drift-Diffusion  Models 

2.1  Mass,  momentum  and  energy  transport  equations 

The  equations  as  presented  here  are  discussed  in  references  [3],  [20],  and  [4].  They  are  derived 
as  the  first  three  moments  of  the  Boltzmann  equation,  with  the  latter  written  for  electrons 
moving  in  an  electric  field  as 

%  +  u-VJ--F-VJ  =  C.  (1) 

at  rn 

Here,  /  =  is  the  numerical  distribution  function  of  a  carrier  species,  x  is  the 

position  vector,  u  is  the  species'  group  velocity  vector,  F  =  F{x^t)  is  the  electric  field,  e 
is  the  electron  charge  modulus,  m  is  the  effective  electron  mass,  and  C  is  the  time  rate  of 
change  of  /  due  to  collisions.  In  the  Boltzmann  equation  above,  it  has  been  assumed  that  the 
traditional  Lorentz  force  field  does  not  have  a  component  induced  by  an  external  magnetic 
field.  The  moment  equations,  which  will  be  derived  subsequently,  are  expressed  in  terms  of 
certain  dependent  variables,  where  n  is  the  electron  concentration,  v  is  the  average  velocity, 
p  is  the  momentum  density,  P  is  the  symmetric  pressure  tensor,  q  is  the  the  heat  flux,  ej  is 
the  internal  energy,  and  Cn,  Cf,,  and  C\v  represent  moments  of  C,  taken  with  respect  to  the 
functions 

hoiu)  =  1, 

hi{u)  =  mil, 


riie  cciuations  ar('  gi'''-'.!  I)y; 


M<0  = 


On 

—  +  V-  (nc)  =  Cn, 


—  -I-  r(  V-  p)  -f  (p  ■  V )(’  =  -(  nF  -V  P  +  Cp. 


(3) 


^  I  i2  .  „  ,  ,mn  ,  ,o  _ 

^(—  1  V  I  +mnei)  +  V-  (v  {—  |  t;  ^  +mne/})  = 

—  env- F  —  V- (vP)  —  V- q  +  Cw-  (4) 

The  first  Maxwell  equation  for  the  electric  potential  must  be  adjoined;  each  species  con¬ 
tributes  a  corresponding  moment  subsystem,  with  appropriately  signed  charge.  We  begin  the 
derivation  with  the  definitions  and  assumptions.  The  concentration  is  given  by  n  :=  /  /  du; 
the  average  velocity  by  u  :=  du;  the  momentum  by  p  :=  mnv;  the  random  velocity 

by  c  :=  u  —  u;  the  pressure  tensor  by  :=  m  f  CiCjf  du;  and  the  internal  energy  density  by 
e/  ^  /  I  c  /  du.  This  function  represents  energy/unit  mass/unit  concentration.  The 
heat  flux  q  is  given  by  :=  y  /cj  |  c  p  /  du.  Finally,  for  reference  in  subsequent  subsec¬ 
tions,  the  electron  current  density  is  given  by  J  :=  —env,  and  the  energy  flux  is  given  by 
S  J  u{y  \  u  \^]f  du.  The  assumptions  on  /  are  now  stated.  The  function  /  is  assumed 
to  decrecise  sufficiently  rapidly  at  infinity: 

lim  hi{u)f{u)  =  0,  f  =  0, 1, 2. 

The  derivation  of  (2),  (3),  (4)  proceeds  by  multiplying  the  Boltzmann  equation  (1)  by 
ho,  hi,  and  /i2,  respectively,  and  integrating  over  group  velocity  space.  With  the  application 
of  certain  standard  identities  ([16]),  the  mass/momentum/energy  system  is  obtained.  In 
addition  to  these  transport  equations,  we  have  Poisson’s  equation  for  the  electric  field,  where 
Ud  :=  doping  and  t  :=  dielectric: 

F  =  -V<^,  (5) 

V-(eV0)  =  -'^ein,-nd.  (6) 

Here,  we  have  used  the  convention  that  there  are  different  species,  each  <  i  concentration  n, 
and  charge  e,.  The  entire  system  consists  of  equations  (2),  (3),  (4)  repeated  according  to 
species,  and  (5),  (6). 

2.2  Moment  closure  and  relaxation  relations 

The  system  derived  in  the  preceding  subsection  has  fifVen  dependent  variables  in  the  case 
of  one  species,  determined  by  (f),  n,  v,  P,  e/,  and  q.  By  moment  closure  is  meant  the 
selection  of  compatible  relations  among  these  vaiiables,  so  that  the  number  of  equations  is 
equal  in  number  to  the  remaining  primitive  variables  selected.  The  relations  to  follow  are 
characterized  by  the  isotropic/parabolic  energy  band  assumption.  We  begin  by  introducing 
a  new  tensor  variable  T,  the  carrier  teuiperaturc,  defined  by 

P,  =  nkl\, 
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where  k  is  Boltzmann’s  constant,  and  a  scalar  variable  W,  the  total  carrier  energy.  A 
program  of  reduction  to  a  set  of  basic  variables,  n,  v,  W ,  and  or  a  set  equivalent  to  these, 
can  be  implemented  by  the  following  assumptions: 


•  The  pressure  tensor  is  isotropic,  with  diagonal  entries  Pg  and  off-diagonal  entries  zero, 
for  a  suitable  scalar  function,  P,.  P,  is  related  to  e;  via  mne/  =  |P,. 

•  It  follows  from  the  previous  assumption  that  temperature  may  be  represented  by  a 
scalar  quantity  T,  and  that  the  internal  energy  is  represented  in  terms  of  T  by 


me/  = 


•  The  total  energy  density  (per  unit  concentration)  w  is  given  by  combining  internal 
energy  and  parabolic  energy  bands  with  m  assumed  constant: 

w  =  mej  -1-  |  v  p, 

and  the  total  energy  (per  unit  volume)  W  is  the  product,  W  =  nw. 


•  The  heat  flux  is  obtained  by  a  differential  expression  involving  the  temperature: 


q  =  —kVT. 

Here,  k  is  the  thermal  conductivity  governed  by  the  Wiedemann-Franz  law  (cf.  [2]), 
described  by 

„  =  +  (7) 

2  e  i  0 

The  standard  choice  for  r  is  r  =  —  1,  but  this  has  some  associated  difficulties.  This 
will  be  amplified  later  in  the  paper.  Here  we  simply  remark  that  the  term  raised  to 
the  exponent  r  in  (7)  is  proportional  to  the  mobility,  which  in  turn  is  proportional  to 
the  momentum  relaxation  time. 

In  the  case  of  N  species,  the  closure  relations  determine  {d  -f  2)N  +  1  variables  in  d 
spatial  dimensions.  It  is  possible  to  rewrite  the  system  (2,  3,  4)  with  the  closure  assumptions 
incorporated.  We  have  the  following. 

=  c„,  (8) 

=  -c7iF~V{nkT)  +  Cp,  (9) 

=  —env-F  —  V-{vnkT) 

+V-{kVT)  +  Cw.  (10) 


dn 

In 


-h  V-  (nv) 


dt 


-f  r(V-  p)  -f  (p  ■  V)n 


dW 

dt 


+  V  -  (  t-  W  ) 
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The  final  step  deals  with  the  replacement  of  the  collision  moments.  Motivated  by  the 
approach  of  [18],  [1],  [20],  and  [11],  we  define  the  recombination  rate  R  and  the  momentum 
and  energy  relaxation  times,  Tp  and  r^,,  respectively,  in  terms  of  averaged  collision  moments 
as  follows. 


1.  The  particle  recombination  rate  R  is  given  by 

R:=-Cn:=-J 

2.  The  momentum  relaxation  time  Tp  is  given  via 


C  du. 


P  [ 

—  :=  —m  I 
Tr,  J 


uC  du  :=  — Cr,. 


3.  The  energy  relaxation  time  is  given  via 

W-Wo  m 


W  -  Wq  rnf  2 

- :=  —  /  I  u  p  /  du  :=  Cw- 

Tyj  2,  J 

Here,  ITo  denotes  the  rest  energy,  |A:To,  where  To  is  the  lattice  temperature. 

The  forms  for  the  relaxation  times  used  in  [1]  and  retained  by  subsequent  authors  are; 


=  Cp(— P, 

Jo 


(11) 

T  +  To  '  2 

Here,  Cp  and  are  physical  constants,  and  the  standard  choice  for  r,  just  as  in  (2.7),  is  —1. 


“  Oil 


T  1 

-  +  -Tp. 


2.3  Drift-diffusion  model 

The  drift-diffusion  model  may  be  obtained  by  taking  zeroth  order  moments  of  the  BTE 
and  adjoining  the  Poisson  equation.  Thus,  one  obtains  the  system  for  N  carriers  with 


concentrations  n,  and  charge  Ci,  i  —  1,  •  ■  • ,  iV: 

^  +  =  -f?.,  (13) 

F  -  -V<^,  (14) 

V-(eVdi)  =  —'^€,n,  —  nd.  (15) 

There  still  remains  the  issue  of  determining  the  constitutive  current  relations.  Classical 

drift-diffusion  theory  gives,  for  N  =  2,  ri)  =  n,  and  n-j  =  p, 

Jn  -  -cpniiVcl)  +  eOnVn,  (16) 

Jp  =  -epppV0  -  eOpVp.  (17) 


0 


The  introduction  of  exponential  relations  for  n  and  p  is  also  common,  as  is  the  use  of  the 
Einstein  relations  linking  the  mobilities,  Pn^Pp,  and  the  diffusion  coefficients  Dn,Dp.  These 


relations  are  specified  by 

r 

Dn  =  {kT/e)fln, 

(18) 

Dp  =  {kTfe)fip. 

(19) 

It  is  also  possible  to  derive  the  constitutive  relations  (16),  (17),  from  the  first  order  moment 
relations  under  the  assumption  that  the  momentum  relaxation  times  tend  to  zero.  The 
details  are  given  in  [20].  In  fact,  the  constitutive  relations  include  a  heat  flux  term  eis  well, 
which  is  suppressed  at  constant  temperature.  If  it  is  not  suppressed,  one  has  an  energy  drift- 
diffusion  model.  In  this  derivation,  one  uses  the  definition  of  mobility  in  terms  of  relaxation 
time. 


3  RT  Models 

In  this  section,  we  shall  employ  a  microscopic  assumption  upon  the  momentum  relaxation 
time,  viz.  ,  we  shall  assume  that  the  collision  term  C  in  (1)  is  of  the  form, 

C  =  (20) 

7p 

where  /i  is  the  odd  part  of  /.  Note  that  this  contrasts  with  the  macroscopic  assumption  on 
Tp,  employed  in  the  hydrodynamic  model  as  described  in  Section  2.  There,  the  representation 
defining  Tp  was  a  post  averaged  expression.  Here,  the  expression  is  employed  in  the  averaging. 
In  this  case,  we  may  obtain  an  expression  for  the  energy  flux  S: 

S  =  -[nfi^F  +  V{nD^)],  (21) 

where  p^  and  are  tensor  expressions  for  mobility  and  diffusion,  defined  in  terms  of 
moments,  and  E  represents  average  energy  per  unit  concentration.  The  details  are  furnished 
in  [5).  It  is  also  shown  there  that  the  current  density  has  the  usual  drift-diffusion  form,  with 
tensor  expressions  for  mobility  and  diffusion.  The  RT  model  makes  the  following  microscopic 
assumptions,  with  distinction  between  £  and  its  average,  E. 

1.  The  even  part  of  /  is  isotropic,  and  a  function  of  £  alone,  and  the  relaxation  time  is 
an  inverse  j)ower  function  of  £\ 


fo  =  fo{£), 

Tp  ==  Tp{£)  =  C£-’-. 
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2.  The  microscopic  kinetic  energy  is  a  quadratic  function  of  £,  and  the  mass  is  not  assumed 


constant: 


G{£) 


where  a  is  an  appropriate  fitting  parameter. 

3.  The  temperature  is  a  modified  variable  in  terms  of  which  the  following  constitutive 
relation  holds  for  E: 

E  =  ^kT{l  +  ^akT).  (23) 

Equation  (23)  allows  for  nonparabolic  energy  bands  as  well  as  non-Maxwellian  distributions. 
Altogether,  the  model  may  be  written  in  terms  of  the  Poisson  equation,  (6),  in  conjunction 
with  the  system, 

V-J  =  0,  (24) 

BF 

V-5  =  J-F-n{—Uu).  (25) 

Here,  J  and  S  have  been  described  previously,  the  latter  in  (21).  In  the  expressions  for  J 
and  5,  the  assumptions  made  for  the  model  lead  to  scalar  representations  for  the  mobility 
and  diffusion  coefficients.  For  example,  the  choice  made  in  [5]  leads  to 

H  =  noTo/T,  (26) 

=  5(1  -  ftrvir.  (27) 

The  diffusion  coefficients  are  defined  by  Einstein’s  relations.  The  collision  term  in  (25)  is 

a  specified  quadratic  function  of  T.  One  significant  advantage  of  the  microscopic  relax¬ 
ation  time  (RT)  assumption  is  that  certain  key  parameters  may  be  fitted  via  Monte-Carlo 
preprocessing,  en  uiring  reliability  of  their  values. 

4  Mathematical  Results  for  the  Hydrodynamic  Model 

In  this  section,  we  shall  describe  some  recent  mathematical  results,  obtained  in  one  spatial 
dimension.  In  the  first  subsection,  we  shall  present  existence  and  boundary  layer  results  for 
a  simplified  version  of  the  steady  state  hydrodynamic  model.  This  will  be  follow'ed  by  a 
convergence  analysis  for  Newton’s  method  in  the  subsonic  case. 


4.1  Existence  and  boundary  layer  theory 

We  first  write  down  the  one  dimensional  evolution  system  in  the  case  of  a  single  carrier,  in 
the  absence  of  recombination. 

dn  ()(nv) 

^  + V  =  0,  (28) 


(29) 


dp  d{pv  -f  knT)  p 

+ - ^ -  =  -enk - , 


dt  dx 

dW  d{vW  +  vknT) 


=  —envF  -f- 


d{K{dT)l{dx))  Vl^  -  Wo 


e-r-  =  -en  -  rid. 
dx 

The  corresponding  steady  state  system  is  obtained  by  setting  the  time  derivatives  equal  to 
zero. 

If  we  rewrite  the  second  steady  state  equation  by  use  of  the  pressure,  P,  we  obtain  the 
equation, 

f )  =  -enF  -  (32) 

dx  Tr, 


The  approach  of  [8]  is  to  eliminate  the  energy  equation  (30)  from  the  system,  and  replace 
its  role  by  a  relation  in  the  spirit  of  gas  dynamics,  i.  e.  by  the  constitutive  relation, 

P{n)  =  KrP,  7  >  1.  (33) 

When  units  are  selected  in  which  e  =  l,e=l,m  =  l,  and  /C  =  1,  we  obtain  the  system,  in 
which  n  and  4)  are  the  only  dependent  variables, 


j  ;=  nv  =  constant, 

(34) 

)x  “ 

(35) 

=  n  ~  nd. 

(36) 

One  can  nominally  specify  boundary  conditions  on  n  and  4>  at  the  endpoints  of  the  device, 
taken  here  as  the  interval,  [0, 1].  If  the  doping  is  such  that  the  built-in  potential  is  the  same 
at  both  ends  of  the  device,  then  we  may  take. 


<^(0)  =  o, 


for  <f>  and 


n(0)  =  no,  n(l)  =  ni,  (38) 

for  n,  where  c^i  must  satisfy  the  following  consistency  relation  with  respect  to  j,  no,  and  ni: 


Here, 


s 


It  is  shown  in  [8]  that  a  weak  solution  exists  for  the  system  (34),  (35),  (36),  satisfying 
the  boundary  conditions  exactly,  or,  in  lieu  of  this,  satisfying  precise  limiting  relationships. 
The  result  for  is  classical,  because  the  equation  is  elliptic.  The  result  for  n  is  provisional, 
and  is  detailed  now.  If  j  >  0,  there  is  a  weak  solution  n  such  that  the  relation, 

n  =  G  +  O',  (41) 

holds,  where  G  is  Holder  (4)  continuous,  and  a  is  monotone  increasing.  Although  n  need 
not  be  of  bounded  variation,  it  is  an  entropy  solution,  in  the  sense  that  the  function  of  a;, 

H{n{x))  -  [F{n{x))  -  F {nmin))signum{n{x)  -  rimin)  +  Cx,  (42) 

is  monotone  increasing.  Here,  Umin  is  a  minimum  (location)  for  F,  and  C  =  sup  S  over 
relevant  arguments.  The  following  precise  statement  is  available  for  subsonic  boundary 
conditions.  If  no,ni  >  Umin,  then  the  following  holds. 

•  Either 

n(l)  =  ni. 


•  or 

lim  n(x) 

X— 

exists,  and  is  a  supersonic  value,  i.  e. ,  is  less  than  rimin,  and  even  less  than  the  conjugate 
value  of  n-i. 

•  Similarly,  either 

n(0)  =  no. 


•  or 


exists  and  is  not  less  than  n„,,„. 


lim  n{x) 

r^O+ 


In  the  second  instance  of  both  cases  above,  boundary  layers  occur,  the  one  on  the  right 
involving  transition  through  the  supersonic  regime. 


4.2  Newton  convergence  theory 

In  this  subsection,  we  order  the  basic  variables  as  v,  n,  T,  and  4>,  because  of  symmetry 
considerations.  Dirichlet  boundary  conditions  are  imposed,  on  n,  T,  4>,  with  n{xm,n)  — 
n(xmax)-  In  one  dimension,  the  steady  state  equations  are  obtained  from  (28),  (29),  (30),  and 
(31)  by  setting  the  time  derivatives  equal  to  zero.  Dirichlet  boundary  conditions  are  imposed 
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in  this  subsection  on  n,  7’,  and  with  n(.r„,m)  =  ni-^max)-  Since,  by  the  conservation  of 
mass  equation,  iiv  =  j,  it  follows  that  the  boundary  conditions  for  v  are  periodic.  The  map 
defined  by  bringing  all  terms  to  *he  left  hand  side  of  the  steady  slate  system  is  called  The 
linearized  eciuations  thus  assume  the  form,  where  the  boundary  conditions  are  homogeneous 
Dirichlet  conditions  for  Sn,  6T,  and  Sq,  and  the  boundary  conditions  on  6v  are  prescribed 
to  be  periodic, 


0 

0 


id- 

n  dr 

—  t 


.  dyr 


) 


+ 


.4  B 
C  D 


6v 

6u 

d 

Sn 

+ 

■  E 

F 

8n 

dx 

6T 

G 

H 

6T 

.  H  . 

(43) 


The  {.spatiall}  )  dependent  eigenvalues  of  the  symmetric  matrix  A  are  calculated  to  be 


Here. 

(45) 

and  the  smaller  eigenvalue  is  positive  if  n  and  T  are  strictly  positive,  and  if 

<  —  ~  c^.  (46) 

m 

This  typ^"  of  [)oint  in  function  si)ace  is  termed  a  sul)sonic  point.  This  case  was  first  considered 
in  [9],  where  damped  .Xewlon/standard  finite  dilTerence  methods  were  presented.  When 
Newton’s  method  is  employed  in  this  way,  it  is  essential  to  determine  conditions  under 
which  the  linear  increments  are  appropriateh  bounded.  This  is  equivalent  to  uniform  bounds 
for  the  operator  derivative  inverse  maps,  and  represents  one  of  the  three  properties  for  an 
(exact)  o]Knator  .Newton  method  to  yield  existence  of  a  root  and  7?-quadratic  convergence. 
The  other  two  are  sullicient  regularity,  and  a  sufficiently  small  starting  residual,  as  measured 
in  tlie  range  space  norm,  b.xplicit  representations  of  B,  C,  D  and  of  F,  G,  H  are  given  in 
[14].  Moreover,  if  the  system  map  <t>,  sulrject  ti  appropriate  Dirichlet  boundary  conditions 
on  11,  r,  and  o.  and  periodic  ijoundary  conditions  on  e,  accordingly  has  the  domain  D<^  C 
A'  =  rii  if  '  ''  X  Hi  ii  arm  range  in  Y  ~  Hi  ^  then  (46)  will  hold  for  every  element 
in  a  closed  ball  /To  C  -V.  centered  at  a  subsonic  element  Uq  £  A'",  such  that  n  and  T  are 
strictly  positive-,  if  is  sufficiently  small.  It  is  appropriate  to  assume  at  the  outset,  then, 
that  D,p  C  Fro  -  every  function  ))oint  in  l),p  satisfies  (46);  we  nray  also  assume  that 

n  and  /'  are  nnilormly  b-ounded  away  Hom  zero  in  this  set. 

rile  l.ipscliitz  jirope-rty  of  the  map  <!>',  wln-re  here  we  view  the  system  (43)  as  the  repre- 
senlatimi  for  T'( 7’.  0)1 -e)  =  /.  with 


A  - 


/? 

V 


V 

r 

Tnn 


;  =  (F/'.,ScT).  /  =  (/,./,). 
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=  O'  r.hii  ). 


is  evident  from  the  representation  for 

The  uniform  inverse  bounding  proceeds  as  follows.  As  shown  in  [9],  with  the  function 
spaces  selected  in  this  paper,  the  //'  product  norm  of  2  can  be  estimated  in  terms  of  the 
norms  of  u.’,  a;',  and  /i,  under  the  conjunction  of  the  hypothesis  (46)  and  the  x 
coerciveness  assumption, 


A  =  /•:  +  E' 


is  uniformly  positive  definite. 


Here,  E"  is  the  matrix  transpose  of  E  and  the  latter  is  defined  by 


an 

+  X  (-At^  • 

.  di  ’  Tp  \mn^  j  dx  . 

A  final  calculation,  making  use  of  the  inner  product  of  [0,u;]  with  (43),  and  the  hypothesis, 

/ill  is  positive,  and  sufficiently  large,  (50) 


where  /in  is  the  nonzero  entry  of  //,  given  by 


,  dv 
hi\  =  - — f- 
dx 


r^-{VV-  VHoX 


with  shows  that  the  lE  product  norm  of  [z,u>]  is  estimated  in  terms  of  the  norm 

of  /.  This  series  of  calculations  controls  the  L°°  norm  of  [2,^];  the  L°°  norm  of  [z',u}"\  is 
now  estimated  by  direct  use  of  the  system  (43),  making  use  of  the  fact  that  [v,n,T,  d>]  € 

We  have  now  outlined  the  proof  of  the  following. 

Theorem  4.1  Let  the  function  spaces  X  and  V  be  selected  as  above,  let  the  steady  state 
system  map  0  be  given  with  Dirichlet  boundary  conditions  on  n,  T ,  4>,  and  periodic  boundary 
conditions  on  v,  such  that  C  where  every  point  in  is  a  subsonic  point,  with  uni¬ 
form  positivity  bounds.  If  (^8)  and  (51)  hold,  and  xq  is  such  that  ‘^(xo)  is  sufficiently  small, 
then  an  R-quadratically  convergent  Newton  sequence  {x„}  may  be  defined  in  the  standard 
way,  with  limit  x,  satisfying  <^(x)  =  0. 


5  Discrete  Schemes  Based  on  Adaptive  Stencils:  ENO 


holds  for  any  real  (  —  (^i,  •  •  •  An  initial  condition  is  adjoined  to  (52). 

For  systems  of  conservation  laws,  local  field  by  field  decomposition  is  used,  to  resolve 
waves  in  different  characteristic  directions.  Analytical  expressions  are  employed  for  the 
eigenvalues  and  eigenvectors  of  an  averaged  Jacobian  matrix.  Typically,  the  Roe  average 
[19]  is  employed.  One  feature  of  the  ENO  schemes  in  [24]  and  [25],  which  is  distinct  from 
the  original  ENO  schemes  of  Harten  et  al  [12],  is  that  multidimensional  regions  are  treated 
dimension  by  dimension:  when  computing  /i(u)x,  in  any  particular  direction,  variables  in  all 
other  directions  are  kept  constant,  and  the  Jacobians  are  treated  in  this  direction.  This,  in 
essence,  reduces  the  determination  of  the  scheme  to  the  case  of  a  single  conservation  law  in 
one  spatial  dimension.  Thus,  to  describe  the  schemes,  consider  the  scalar  one  dimensional 
problem,  and  a  conservative  approximation  of  the  spatial  operator  given  by 


Here,  the  numerical  flux  /  is  assumed  consistent: 


/j  +  L  =  /(«;-(,•••  ,u_,+fc);  /(u,---,u)  = /(u).  (54) 

The  conservative  scheme  (53),  which  characterizes  the  /  divided  difference  as  an  approx¬ 
imation  to  suggests  that  /  can  be  identified  with  an  appropriate  function  h  satisfying 

f{u{x))  =  r  J  MO  dO  (55) 

If  //  is  any  j)!imitive  of  h,  then  h  can  be  computed  from  //'.  H  itself  can  be  constructed 
hy  .Newton’s  di\  ided  difference  method,  beginning  with  differences  of  order  one.  since  the 
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constant  term  is  arbitrary.  The  necessary  divided  differences  of  H,  of  a  given  order,  are 
expressed  as  constant  multiples  of  those  of  /  of  order  one  lower.  After  the  polynomial  Q  of 
degree  r  +  1  has  been  constructed,  set 


(56) 


to  obtain  an  rth  order  method.  The  construction  is  based  on  an  adaptive  stencil  in  the 
following  sense; 


•  One  begins  with  an  appropriate  starting  point  to  the  left  or  right  of  the  current  “cell” 
by  means  of  upwinding  as  determined  by  the  sign  of  the  derivative  of  a  selected  flux. 

•  As  the  order  of  the  divided  differences  is  increased,  the  divided  differences  themselves 
determine  the  stencil:  the  “smaller”  divided  difference  is  chosen  from  two  possible 
choices  at  each  stage,  ensuring  a  smoothest  fit. 

•  Lax-Friedrichs  building  blocks  or  Roe  building  blocks  can  both  be  used.  For  the  lat¬ 
ter,  in  cells  with  sonic  points,  a  local  Lax-Friedrichs  building  block  is  used  to  avoid 
expansion  shocks. 


Steady  states  are  reached  by  explicit  time  stepping  of  arbitrary  order;  nonstandard  high 
order  Runge-Kutta  methods  exist  [24]  which  preserve  nonlinear  stability  of  the  first  order 
Euler  forward  version  under  suitable  CFL  time  step  restrictions.  The  computer  program 
is  fully  vectorized  for  computations  on  Cray  supercomputers.  For  details  of  the  efficient 
implementation,  see  [23]. 


6  Conservation  Law  Format  for  Hydrodynamic  and 
RT  Models 

In  this  section,  we  shall  specify  the  conservation  law  format  for  the  two  dimensional  hydro- 
dynamic  model,  and  for  the  one  dimensional  RT  model. 

6.1  Hydrodynamic  model  conservation  format 

Define  the  vector  of  dependent  variables  as 

u  =  {n,a,T,W),  (57) 

where  p  =  (cr, r).  The  system  (8),  (9),  (10)  can  be  written  in  the  concise  form,  in  two 
dimensions,  as 

ut  +  /i(«0x  +  /2(u)y  =  c(«)  +  G(u,<j!>)  +  (0,0,0,  V-  (kVT)).  (58) 
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The  following  identifications  have  been  made  in  (58). 


ft  ^ 

^  m’  3  mn  2mn  ’  mn’  3mn  3m^n^  ’ 

T  err  2  cr^  5tW  cr^  + 

^  m’  mn’  3  mn  2mn  ’  3mn  3m^n'^  ’ 


(59) 

(60) 


c(u) 

G(u) 


W-Wq 

Tw 


), 


(0,  —e7iFi,  — en/2,  -enF-  v). 


(61) 

(62) 


The  eigenvalues  and  eigenvectors  of  f[  and  are  known  (cf.  [23]),  and  are  readily 
incorporated  into  the  field  by  field  decomposition  required  for  the  implementation  of  ENO. 


6.2  RT  conservation  format 

We  shall  present  the  conservation  law  form  of  the  RT  model.  We  begin  with  the  vector  form, 

ut  +  f{u)x  =  9{u)xx  + Hu).  (63) 


In  equation  (63), 

nE 

u  =  (en, - ), 

m 

f{u)  =  4>'n{eHE),  ti^{E)  +  D{E)), 
g{u)  =  (nD(E),  nD^{E)), 

p\  rp 

Hu)  =  {0,  enHE){<i>')'^  +  ^{n  -  nd)nD{E)  -  n{—  \coii))- 


(64) 

(65) 

(66) 
(67) 


It  can  be  shown  that  the  left  hand  side  defines  a  hyperbolic  system,  since  the  eigenvalues  of 
f'{u)  are  real,  for  all  positive  n  and  T. 


7  Numerical  Simulation  Results 

We  now  present  numerical  simulation  results  for  one  carrier,  two  dimensional  MESFET 
devices  and  one  dimensional  diodes.  The  third  order  ENO  shock-capturing  algorithm  with 
Lax- Friedrichs  building  blocks,  as  described  briefly  in  Section  5  and  in  more  detail  in  [25],  is 
applied  to  the  hyperbolic  part  (the  left  hand  side)  of  Equations  (6.2)  and  (6.7).  A  nonlinearly 
stable  third  order  Ilunge-Kutta  time  discretization  [24]  is  used  for  the  time  evolution  towards 
steady  states.  The  forcing  terms  on  the  right  hand  side  of  (6.2)  and  (6.7)  are  treated  in  a 
time  consistent  way  in  tlu;  Runge-Kutta  time  stepping.  The  double  derivative  terms  on  the 


right  hand  side  of  (6.2)  and  (6.7)  are  approximated  by  standard  central  differences  owing  to 
their  dissipative  nature.  The  Poisson  equation  (2.6)  is  solved  by  direct  Gauss  elimination  for 
one  spatial  dimension  and  by  Successive  Over-Relaxation  (SOR)  or  the  Conjugate  Gradient 
(CG)  method  for  two  spatial  dimensions.  Initial  conditions  are  chosen  a.s  n  —  rid  for  the 
concentration,  T  =  To  for  the  temperature,  and  u  —  v  =  0  (two  spatial  dimensions)  or  u  =  0 
(one  spatial  dimension)  for  the  velocities.  A  continuation  method  is  used  to  reach  the  steady 
state:  the  voltage  bias  is  taken  initially  as  zero  and  is  gradually  increased  to  the  required 
value,  with  the  steady  state  solution  of  a  lower  biased  case  used  as  the  initial  condition  for 
a  higher  one. 

7.1  Two  dimensional  MESFET 

We  simulate,  using  the  Hydrodynamic  model  (6.2)-(2.6),  a  two  dimensional  MESFET  of  the 
size  0.6  X  0.2/rm^.  The  source  and  the  drain  each  occupies  0.1/rm  at  the  upper  left  and  the 
upper  right,  respectively,  with  the  gate  occupying  0.2/im  at  the  upper  middle  (Figure  1,  left). 
The  doping  is  defined  by  uj  =  3  x  in  [0,0.1]  x  [0.15, 0.2]  and  in  [0.5,  0.6]  x  [0.15, 0.2], 

and  Hd  =  1  X  10''c?7i“^  elsewhere,  with  abrupt  junctions  (Figure  1,  right).  A  uniform  grid 
of  96  X  32  points  is  used.  Notice  that  even  if  we  may  not  have  shocks  in  the  solution,  the 
initial  condition  n  =  rid  is  discontinuous,  and  the  final  steady  state  solution  has  a  sharp 
transition  around  the  junction.  With  the  relatively  coarse  grid  we  use,  the  noii-oscillatory 
shock  capturing  feature  of  the  ENO  algorithm  is  essential  for  the  stability  of  the  numerical 
procedure. 


Figure  1:  Two  dimensional  MESFET.  Left;  the  geometry;  Right;  the  doping  /jj 
We  apply,  at  the  source  and  drain,  a  voltage  bias  vbias  =  2V^.  The  gate  is  a  Schottky 
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contact,  with  a  negative  voltage  bias  vgnte  =  — U.St'’  and  a  V(>ry  low  ((jiufiilration  value 
n  =  3.9  X  10® cm~^  obtained  from  Equation  (5.1-19)  of  [23],  Tl’.c  lattice  temperature 
taken  as  Tq  =  300° /t'.  The  numerical  boundary  conditions  are  summarized  as  follows  (where 
$0  =  ^  In  with  kb  =  0.138  x  10'“*,  c  =  0.1602,  and  n,  =  1.4  x  in  our  units); 

•  At  the  source  (0  <  x  <  0.1,y  =  0.2);  <1>  =  Oq  for  the  potential;  n  =  3  x  lO^'cm”^ 
for  the  concentration;  T  =  300°  A'  for  the  temperature;  u  =  Ofnnjps  for  the  horizontal 
velocity;  and  Neumann  boundary  condition  for  the  vertical  velocity  v  (i.e.  ||  =  0 
where  n  is  the  normal  direction  of  the  boundary). 

•  At  the  drain  (0.5  <  x  <  0.6,  y  =  0.2):  $  =  ^>o  +  vbias  =  +  2  for  the  potential; 

n  =  3  X  10^ 'cm~^  for  the  concentration;  T  =  300° K  for  the  temperature;  u  =  Opm/ps 
for  the  horizontal  velocity;  and  Neumann  boundary  condition  for  the  vertical  velocity 

V. 

•  At  the  gate  (0.2  <  x  <  0,4,  y  =  0.2):  =  <&o  +  vgate  =  ^>o  —  0.8  for  the  potential; 

n  =  3.9  X  10®cm“^  for  the  concentration;  T  =  300° K  for  the  temperature;  u  =  0pm fps 
for  the  horizontal  velocity;  and  Neumann  boundary  condition  for  the  vertical  velocity 

V. 

•  At  all  oilier  parts  of  the  boundary  (0.1  <  x  <  0.2,  y  =  0.2;  0.4  <  x  <  0.5,  y  =  0.2; 
X  =  0,0  <  y  <  0.2;  x  =  0.6,0  <  y  <  0.2;  and  0  <  x  <  0.6, y  =  0),  all  variables  are 
equipped  with  Newmann  boundary  conditions. 

The  boundary  conditions  chosen  are  based  upon  physical  and  numerical  considerations. 
They  may  not  be  adequate  mathematically,  as  is  evident  from  some  serious  boundary  layers 
observable  in  Figures  2  through  6.  ENO  methods,  owing  to  their  upwind  nature,  are  robust 
to  different  boundary  conditions  (including  over-specified  boundary  conditions)  and  do  not 
exhibit  numerical  difficulties  in  the  presence  of  such  boundary  layers,  even  with  the  extremely 
low  concentration  prescribed  at  the  gate  (around  10~'^  relative  to  the  high  doping).  We  point 
out,  however,  that  boundary  conditions  affect  the  global  solution  significantly.  We  have  also 
simulated  the  same  problem  with  different  boundary  conditions,  for  example  with  Dirichlct 
boundary  conditions  everywhere  for  the  temperature,  or  with  Neumann  boundary  conditions 
for  all  variables  except  for  the  potential  at  the  contacts.  The  numerical  results  (not  shown 
in  this  paper)  are  noticeably  different.  This  indicates  the  importance  of  studying  adequate 
boundary  conditions,  bom  both  a  physical  and  a  mathematical  point  of  view. 

In  Figures  2  through  6,  we  show  pictures  of  the  concentration  n,  temperature  T,  horizontal 
velocity  u,  vertical  velocity  e,  and  the  potential  <I>.  Surfaces  of  the  solution  are  shown  at 
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the  left,  and  cuts  at  j/  —  0.175,  which  cut  through  the  middle  of  the  high  doping  “blobs” 
horizontally,  are  shown  at  the  right. 


Figure  2:  Two  dimensional  MESFET,  concentration  n.  Left:  surface  of  the  solution; 
Right:  cut  hi  y  ~  0.175 


Figure  3:  Two  dimensional  MESFET,  temperature  T.  Left:  surface  of  the  solution; 
Right:  cut  hi  y  =  0.175 
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Figui 
Right:  c 


Figu 


Figure  6:  Two  dimensional  MESFET,  potential  Left;  surface  of  the  solution;  Right: 
cut  at  ?/  =  0. 175 

Notice  that  there  is  a  boundary  layer  for  the  concentration  n  at  the  drain  but  not  at  the 
source.  Also  notice  the  rapid  drop  of  n  at  the  depiction  region  near  the  gate.  The  temperature 
achieves  its  maximum  around  the  left  corner  of  the  drain.  The  leakage  current  at  the  gate 
appears  negligible  from  the  normal  velocity  component,  while  the  horizontal  component 
shows  evidence  of  strong  carrier  movement  toward  the  source  beneath  the  left  gate  area,  and 
strong  movement  toward  the  drain  immediately  to  the  left  of  the  drain  junction. 

We  have  also  simulated  the  same  MESFET  with  a  higher  doping  ratio:  3  x  in 

the  high  doping  region  versus  1  x  10*'^c/n“^  in  the  low  doping  region.  We  observe  similar 
results  (pictures  not  shown  here). 

7.2  HD  model  for  a  one  dimensional  diode  —  spurious  velocity 
overshoot 

A  notorious  phenomenon  of  HD  models  is  that  spurious  velocity  overshoot  occurs  at  the 
drain  junction  of  an  diode.  It  is  intrinsic  to  the  model  and  is  not  a  numerical 

artifact,  as  is  verified  by  our  grid  refinement  study  and  by  using  different  numerical  algo¬ 
rithms.  This  phenomenon  is  closely  related  to  the  physical  assumption  governing  the  heat 
conduction  term.  (Iiiudi.  Odeh  and  Rudan  [10]  observed  that  the  spurious  overshoot  can  be 
greatly  reduced  by  an  empirical  modification  of  the  Wiedemann-Franz  law  for  the  thermal 
conductivity. 

In  this  subsection  we  perform  an  extensive  numerical  study  of  the  dependency  of  the 
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spurious  velocity  ovevslioot  upon  the  heat  conduction  term.  The  diode  we  simulate 

has  a  length  0.6/^//^  with  a  doping  defined  by  nj  =  3x  in  [0, 0.1]  and  in  [0.5, 0.6],  and 

Tij  =  1  X  10‘®cni“^  in  [0.15,0.45],  with  smooth  junctions  (Figure  7).  The  lattice  temperature 
is  taken  as  To  —  29G.‘21°/\  .  We  apply  a  voltage  vbias  =  1.5V'',  as  is  the  case  in  [10]. 
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Figure  7;  The  T nj 


for  the  one  dimensional  n'^-n-n'^  diode 


Figure  8:  HD  for  one  dimensional  n^-n-n'^  diode.  Velocity  u  (upper),  temperature  T 
(middle  left),  concentration  n  (middle  right),  total  energy  W  (lower  left)  and  electric  field 
— (lower  right).  Solid  line:  r=-l;  dashed  line:  r=-2;  circles:  r=-2  in  the  coefficient  of  (2.7); 
pluses:  r=-2  in  (2.7) 
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The  standard  IlD  model  uses  r  —  —1  in  the  Wiedemann- Franz  law  (2.7)  and  the  relax¬ 
ation  times  (2.11).  The  numerical  solution  of  this  is  shown  as  solid  lines  in  Figure  8.  We 
can  clearly  observe  the  spurious  velocity  overshoot  at  the  right  junction,  but  otherwise  the 
solution  is  basically  correct  comparing  with  direct  Monto-Carlo  simulations  (not  shown). 
When  r  is  taken  as  —2  in  (2.7)  and  (2.11),  the  solution  is  completely  wrong  (dashed  line  in 
Figures  8).  However,  when  one  takes  r  =  —2  only  in  the  coefficient  of  k  in  (2.7)  but  leaves 
r  —  —1  in  the  power  of  k  in  (2.7)  and  in  (2.11),  i.e.,  when  one  uses 

„  =  (68) 
2  e  Iq 

in  the  place  of  (2.7)  and  leaves  r  =  —1  in  (2.1 1)  unchanged,  as  was  done  in  [10],  one  obtains  a 
greatly  reduced  spurious  overshoot  (the  circles  in  Figures  8).  Finally,  the  result  with  r  =  —2 
in  K  in  (2.7)  but  with  r  =  —  1  in  (2.11)  unchanged,  is  shown  by  pluses  in  Figure  8.  We  can 
see  that  the  spurious  overshoot  also  disappears. 

7.3  RT  model  for  a  one  dimensional  diode 

We  present  numerical  simulation  results  for  the  RT  model,  described  in  Section  3,  for  the 
same  one  dimensional  diode  used  in  Subsection  7.2.  Although  the  RT  model  is  a  parabolic 
system  with  two  equations,  the  existence  of  sharp  transition  regions  near  the  junctions 
justifies  the  usage  of  FNO  shock  caj)turing  algorithms  for  the  hyperbolic  part. 

In  Figure  9,  we  show  the  results  of  velocity  u,  temperature  T,  concentration  n,  total 
energy  E,  and  electric  lic'ld  of  the  RT  simulation,  in  circles,  in  a  background  of  standard 
HD  results  {r  -  -1)  in  solid  lines,  and  of  HD  results  with  r  =  -2  in  the  coefficient  of  k  in 
(2.7)  but  with  r  -  -1  in  tlu'  power  (i.e.,  (2.7)  is  replaced  by  (7.1)),  and  r  =  — 1  in  (2.11),  in 
dashed  line.^.  W'r  can  see  that  the  RT  model  greatly  reduces  the  spurious  velocity  overshoot 
and  is  compiUable  with  the  result  of  the  empirically  modified  HD  result  in  dashed  lines. 

Exten.sive  numerical  tests  about  the  RT  model,  as  well  as  comparisons  between  the  RT 
and  HD  models,  (  institute  ongoing  research,  jointly  with  U.  Ravaioli,  E.  Kan  and  D.  Chen 
at  ' ;  '  d  w  h'  [ihi.ois. 
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Figure  9;  One  dimensional  n^-n-ri^  diode.  Velocity  u  (upper),  temperature  T  (middle 
left),  concentration  n  (middle  right),  total  energy  E  (lower  left),  and  electric  field  (lower 
right).  Solid  line;  standard  HD  with  r=-l;  circles;  RT;  dashed  lines:  HD  with  r=-2  in  the 
coefficient  of  (2.7) 
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